The Higgs mode in a two-dimensional superfluid 
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We present solid evidence for the existence of a weh-defined Higgs amphtude mode in two- 
dimensional relativistic field theories based on analytically continued results from quantum Monte 
Carlo simulations of the Bose-Hubbard model in the vicinity of the superfluid-Mott insulator quan- 
tum critical point, featuring emergent particle- hole symmetry and Lorentz-invariance. The Higgs 
boson, seen as a well-defined low-frequency resonance in the spectral density, is quickly pushed to 
high energies in the superfluid phase and disappears by merging with the broad secondary peak 
at the characteristic interaction scale. Simulations of a trapped system of ultra-cold ^^Rb atoms 
demonstrate that the low-frequency resonance is lost for typical experimental parameters, while the 
characteristic frequency for the onset of strong response is preserved. 
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The emergence of low-energy excitations in systems 
with spontaneously broken symmetry is one of the most 
fascinating and fundamental subjects in physics relevant 
for understanding such diverse systems as solids, mag- 
nets, ultra-cold atoms, and relativistic fields. The gener- 
ation of mass by the Anderson-Higgs mechanism (TJ [2] is 
particularly important for the Standard model [3], where 
detection of the Higgs boson is still the missing link in 
revealing this mechanism, as well as for numerous super- 
fiuid/superconducting condensed- matter systems. In re- 
alistic materials the amplitude mode is often masked by 
other low-energy excitations. These complications are 
avoided by considering atomic bosonic superfluids, de- 
scribed by a complex order parameter field, which con- 
stitute the cleanest experimental realization. 

Generic superfluids do not feature a well-defined Higgs 
boson (by 'well-defined' we understand a mode seen as a 
sharp resonance). Weakly-interacting gases do not have 
it because at and around the critical temperature Tc for 
the superfluid (SF) to normal fluid phase transition all 
long-wave elementary excitations are overdamped, while 
at T ^ the low-energy spectrum is exhausted by the 
Bogoliubov quasiparticle excitations where phase and 
density are canonical variables. Strong interactions do 
not necessarily change this picture. As long as the critical 
temperature remains large, as in ^He, long-wave excita- 
tions are overdamped at \T — Tc\ <^Tc. At low tempera- 
ture, only the Nambu-Goldstone phase modes remain at 
low frequencies whereas the amplitude mode is pushed 
to the incoherent continuum at the (large) characteris- 
tic interaction energy scale. Suppressing Tc by increas- 
ing interactions may trigger a first order transition to 
the solid phase and not work either. It is thus crucial 
to consider an experimental system with a second-order 
quantum critical point (QCP) where Tc for superfluidity 
can be tuned to near zero. 



The Bose-Hubbard model 
H = -JY1 blb^ + ^Y.''^{n,-^)-J2{f^-Vi)ni, (1) 

<ij> i i 

with experimentally adjustable ratios between the hop- 
ping amplitude J, on-site interaction [/, chemical po- 
tential /i, and trapping potential v^, provides an accurate 
description of ultra-cold bosonic atoms in optical lattices. 
At integer filling factor, = (n^), and zero temperature it 
undergoes a second-order quantum phase transition from 
SF to the Mott insulator (MI) phase as the interaction 
strength is increased [4]. The critical field theory be- 
hind this transition is Lorentz- invariant and particle-hole 
symmetric (while the SF-MI transition for generic values 
of u belongs to the universality class of the ideal Bose 
gas at vanishing density and is excluded) |4| i5j . Despite 
the decay into two phase modes the existence of a sharp 
Higgs boson is guaranteed in two special limits: (i) in a 
three-dimensional (3D) system where the corresponding 
4D quantum field theory is at the upper critical dimen- 
sion with asymptotically exact mean-field behavior and 
vanishing decay rates (see Ref. [6 for an experimental 
observation in a quantum antiferromagnet); (ii) at large 
momentum when the relativistic time dilation effect leads 
to an increased quasiparticle decay time. 

The most intriguing question is whether the low- 
frequency Higgs boson can be seen as a well-defined ex- 
citation at zero momentum at the density-driven QCP of 
the 2D Bose-Hubbard model and how it disappears with 
detuning to the SF phase. Equally important are finite 
temperature effects and the role of the trapping poten- 
tial in experiments. A theoretical treatment of the Higgs 
amplitude mode is notoriously difficult and controversial. 
In Refs. [51 [THS] exact scaling laws in the low-frequency 
limit were established, as well as arguments given that 
the mode is at the edge of the two-phonon continuum, 
rendering the mode overdamped. Huber et at. used a 
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FIG. 1. Universal scaling predictions for the scalar suscepti- 
bility (solid lines). The dashed-dotted line depicts prediction 
of Ref. [12^ and misses most of the spectral density at the 
relevant energy scale A oc (1 — U /UcY . The two alternatives 
for connecting universal power laws are shown by dashed lines 
(one may also imagine multiple peaks in the crossover region) . 



variational Ansatz which, however, predicted a spurious 
first order SF-MI transition, and thus was limited to the 
parameter regime away from quantum criticality ^3 ^ . 
Podolsky et al. generalized the field theoretical results of 
Ref. [7] to high frequencies and discussed in detail the re- 
sponse function for the order-parameter density (scalar 
response) within 1/N and a weak coupling expansion 
schemes [12]. They revealed a broad peak whose maxi- 
mum saturates at finite value at the QCP and concluded 
that close enough to the transition, it becomes impossible 
to identify the Higgs energy with the peak position. Their 
findings are in quantitative and qualitative disagreement 
with those reported here, as is detailed in the supplemen- 
tary material [13]. The major problem with the results 
of Ref. [12] is the strong violation of the universal low- 
frequency scaling law for the scalar response function [5^ , 
S{uj) (X A3-2/^F(cj/A), where A ex (1 - U/UcY is the 
characteristic energy scale in proximity of the quantum 
critical point, and v = 0.6717 the correlation length ex- 
ponent. As is shown in Fig.[l] the theory misses most of 
the spectral density in the A < a; < 4 J range. 

In this Letter, we employ quantum Monte Carlo sim- 
ulations of the 2D model ([T]) in the lattice path integral 
representation using the worm algorithm [TIHTB] to study 
the spectral density of the kinetic energy correlation func- 
tion at zero momentum, in combination with an analytic 
continuation method. We unambiguously demonstrate 
the existence of a low-energy resonance peak associated 
with the Higgs boson in close vicinity of the QCP by dis- 
criminating it from the second broad peak at the typical 
lattice-model energies. The Higgs boson energy, cjr, ob- 
tained from the peak maximum increases with detuning 
nearly identically as that of the particle-hole gap Ami in 
the MI phase. The spectral density associated with the 



Higgs boson broadens with detuning and quickly over- 
laps with other higher energy modes: It is no longer seen 
as a resonance peak for a detuning as small as 20 %, in 
line with the parameter regime where particle and hole 
masses were found to be equal on the MI side [17 . On the 
other hand, in close vicinity of the QCP the Higgs boson 
remains visible in the spectral density at temperatures 
as high as Tc. A peak is even seen in the normal phase; 
only at a temperature T > 2Tc the Higgs resonance is no 
longer visible. However, the onset of strong response at 
low-frequency is barely modified. These results, further 
supported by simulations of realistic trapped systems, 
explain why the experimental protocol of extracting cjh 
from the onset of strong response [18^ works even in the 
absence of low-frequency resonance. 
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FIG. 2. Spectral density S{(jj) of the kinetic energy correla- 
tion function for U / J — 16 (thick solid line), 14 (dashed line), 
and 12 (thin solid line) at low temperature T/J — 0.1. The 
Higgs amplitude mode (cjh) emerges as a well-defined peak 
on approach to the quantum critical point at Uc = 16.7424. 

A small uniform modulation of the optical lattice 
depth [19 leads, under mapping to the Bose-Hubbard 
model Q, to a perturbation proportional to the total 
kinetic energy of the system [20j, K 



V = S{t)K, S{t) = 



SJjt) 
J 



(2) 



where the small SJ{t) is proportional to the lattice mod- 
ulation amplitude. Within standard linear response the- 
ory one computes the corresponding correlation function 
x{ioJn) = {K{r)K{0))i^^ + (K) at Matsubara frequen- 
cies, ujn = 27rTn, and performs an analytic continuation 
to obtain its spectral density S{uj). This quantity is di- 
rectly proportional to the energy absorbed by the system 
in the experiment [TF . In the path integral represen- 
tation K{iujn) has a straightforward Monte Carlo esti- 
mator, ^j.e'^^'^^^^ where the sum goes over all hopping 
transitions in the imaginary time evolution of the system. 
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Thus x(ico'n) is computed exactly, i.e., the error bars are 
statistical and they can be reduced arbitrarily by increas- 
ing the simulation time. Our relative error bars for the 
lowest frequencies are of the order of 10~^. Nevertheless, 
long simulations are required because finite x is found 
only after the cancelation of macroscopic factors. The 
combination of the linear system size (L = 20 in prac- 
tice) of our square lattice, temperature, and U/T has to 
be chosen such that L is always significantly larger than 
the correlation length by a factor of at least four in order 
for finite size effects to be negligible. In the supplemen- 
tary material we provide details about the analytic con- 
tinuation procedure which extracts the spectral density 
S{uj) = Imx(cc;) from x(^^n) and show tests confirming 
the reliability of the results reported here [13 . 

Unambiguous evidence for the existence of the sharp 
amplitude mode in the vicinity of the quantum critical 
point located at {U/J)c = 16.7424 is provided in Fig. 
[21 We observe a well-defined resonance in the normalized 
spectral density at a scale much smaller than [/, which 
softens on approach to the QCP in a way that is com- 
patible with the 3D XY universality class, as is shown 
in Fig. [3) We identify the energy of the Higgs boson 
with the peak maximum. The width of the resonance 
peak narrows when U ^ Uc, suggesting that this peak 
is part of the universal scaling scaling function [5 which 
can be rewritten as S{(jj) oc cj^ '^^^ F{ijO /uoh)- Since at 
frequencies uoh ^ ^ the response must be independent 
of ujh we have S{uj > uh) oc a;^-^/^ = a;^-0225^ ^ 
is increasing extremely slowly. The overall picture in the 
asymptotic limit is that of a Higgs peak superimposed 
on a smeared step function, see Fig. [l] When tuning 
away from the critical point, the Higgs mode broadens 
and overlaps with the second peak around the crossover 
scale U/J^ 12. Beyond this point the Higgs boson can 
no longer be discerned as a separate mode, as is shown 
in Fig. |2j This limits the observation of the amplitude 
mode to the region in close vicinity of the QCP. 

We also would like stress that the resonance at ujh is 
seen on both sides of the QCP, i.e. it is also present in the 
MI phase close to the QCP. The similarity between the 
two responses for the same amount of detuning from the 
critical point, as is evident from Fig.[4j is expected inside 
the correlation volume despite obvious differences in the 
low-frequency part. This result provides further evidence 
that the analytic continuation procedure is stable. 

Though our imaginary time data decaying as ~ 1/r^ 
for large r is compatible with the scaling prediction [71 [8] 
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for 



0, our errors bars at large imaginary 



times are too large and our system sizes too small to 
resolve this law unambiguously in analytic continuation. 

There is substantial room for incoherent spectral 
weight between cjr and u ^ U which can be filled by 
other modes such as 'doublon' (double occupancy) ex- 
citations, pairs of phase modes with zero total momen- 
tum, as well as multi-Higgs modes. Our interpretation 
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FIG. 3. (Color online). Characteristic energies in the vicin- 
ity of the quantum critical point at {U/J)c = 16.7424. Black 
circles for U > Uc and U < Uc stand for particle-hole gaps 
in the MI phase [17] and energies of the Higgs bosons, re- 
spectively. Red squares denote the location of the broad 
secondary peak in S{uj) until it merges with the amplitude 
mode at interaction strength U < 12J to form a single peak 
(blue squares). Shaded regions indicate the characteristic 
broadening of peaks. The thick black line is the critical law 
2.25 J\{U — Uc)/ obtained by fitting the smallest MI gaps; 
its mirror reflection is shown as thin black line. The dashed 
red line indicates the typical interaction scale U/ J. 



S(C0)| ASF [U/J=16.5] 




CO /J 



FIG. 4. Spectral densities in the SF and MI phases approx- 
imately with the same amount of detuning from quantum 
criticality \Uc - U\/U ^ 0.015. 



of the data is that higher frequency doublons ('screened' 
by interaction effects) overlap with lower frequency crit- 
ical phase modes creating an intermediate broad peak at 
frequencies between U and z J, except extremely close to 
QCP where the second peak saturates at about 6 J when 
tuning U Uc We associate it with pairs of phase 
modes with opposite momenta near the Brillouin zone 
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boundary (which dominate in the integral over momen- 
tum space). This classification, however, is not rigorous 
in the quantum critical region [5 , as is evidenced by the 
similarity between the SF and MI responses. 

Current experiments with ultra-cold atoms are typi- 
cally performed at finite temperature T /U > 0.05 (such 
that T/Jr^ C^(l) at = Uc) and in the presence of 
parabolic confinement. In Fig. [5] we demonstrate that 
for the representative case U/J = 16 with Tc/J ~ 0.45 
the Higgs mode remains clearly visible at all tempera- 
tures below the superfiuid transition temperature and 
even slightly above it! At a temperature T > 2Tc the 
two peaks finally merge . Nevertheless, 5'(cj), still levels 
off at the amplitude mode frequency uo ^ ujyl and has the 
same frequency for the onset of strong response, meaning 
that these features can be used to extract ujh experimen- 
tally. Note that phase coherence can extend across finite 
systems at temperatures well above the thermodynamic 
Tc for the Kosterlitz-Thouless transition characterized by 
an exponentially divergent correlation length. 



U/J=16 

T/J=0.1, 0.2, 0.5, 1 




FIG. 5. (Color online). Evolution of spectral density with 
temperature at U/J= 16. As temperature is increased from 
T/J = 0.1 (thin black line), to T/J = 0.2 (thin dashed line), 
and T/J = 0.5 (thick blue line), the peaks get broader but 
remain clearly identifiable. At T/J = 1 (thick dashed green 
line) the two peaks merge. 

Inhomogeneous broadening of spectral density caused 
by the trap has a dramatic effect on the structure of 
S{uj) as signals from different parts of the system are su- 
perimposed on each other. Moreover, in the presence 
of external potential gradients the spectral density is 
no longer vanishing at ^ because of low-frequency 
sound modes (predominantly in the trap edges), in line 
with experimental observations JQ. Under these condi- 
tions, the Higgs mode can no longer be seen as a sharp 
resonance in S{uj). There is instead a broad maximum 
with irregular shape. Finite temperature effects further 
transform it into a smooth single peak. In Fig.[6]we show 



the comparison between the homogeneous and trapped 
cases. The simulation was performed for realistic exper- 
imental parameters [18] but at a variety of temperatures 
in order to discriminate between trap and temperature 
effects. Even though the resonance is lost in the total sig- 
nal, the steep onset of spectral response still correlates 
remarkably well with the energy of the Higgs boson. 



S(co) 



T/J=0.5 




FIG. 6. Effects of trapping potential and finite temperature 
for U/J = 14 (Tc/J ^ 1.04). The spectral density of a ho- 
mogeneous system at low temperature with Higgs resonance 
(dashed line) transforms into a broad (irregular) peak due to 
inhomogeneous broadening in a trapped system with N — 190 
particles at unit filling factor in the middle (thin solid line). 
At a temperature Tj J — 0.5 (thick solid line) we observe a 
smooth single peak. 



In conclusion, we find that the Higgs boson is a well- 
defined though significantly damped excitation in close 
vicinity (~ 20%) of the particle-hole symmetric and 
Lorentz- invariant SF-MI quantum critical point in two 
dimensions. It is seen as a resonance in the kinetic energy 
correlation function which is directly probed through 
the modulation of the optical lattice depth in experi- 
ments with ultra-cold atoms. The energies of the am- 
plitude mode match particle-hole gaps in the Mott in- 
sulator phase for the same amount of detuning away 
from quantum criticality. While temperatures at least 
as high as the critical temperature for superfiuidity pre- 
serve the Higgs resonance, inhomogeneous broadening in 
small trapped systems erases resonance-type features in 
the spectral function. Nevertheless, it is possible to de- 
termine the energy of the amplitude mode from the onset 
of strong response, as is done in a recent experiment J8]. 

This work was supported by the National Science 
Foundation under grant PHY- 1005 543. Use was made 
of the maxent application [22] in the ALPS libraries [23] . 
We thank S. Sachdev, B. Svistunov, W. Zwerger, and the 
authors of Ref. [18 for insightful discussions, as well as 
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FIG. 7. (Color online). Kinetic energy autocorrelation func- 
tion for L = 20, T/J = 0.1, and U/J = 16.5 (black circles) 
described by the near exponential decay with integrated au- 
tocorrelation time ta ~ 0.5 sec. 

for sharing their results with us prior to publication. 



Supplementary material for 
"The Higgs mode in a two- 
dimensional superfluid" 

QUANTUM MONTE CARLO SIMULATIONS 

We simulate the two-dimensional Bose-Hubbard model 

<ij> i i 

with hopping amplitude J, on-site interaction [/, chemi- 
cal potential /i, and trapping potential Vi using quantum 
Monte Carlo simulations in the path-integral representa- 
tion with worm-type updates p^HTB] . Unless otherwise 
stated, our unit is the hopping J = 1. The worm algo- 
rithm is a well established method in the study of bosonic 
systems such as cold gases and Helium-4. It can not only 
simulate bosons loaded in an optical lattice realistically, 
but even go beyond that and study a million particles 
for temperatures of the order of the hopping on a single 
CPU. It is well established that Worm Algorithm does 
not suffer from the critical slowing down problem and 
remains efficient on approach to the critical point. For 
our system size at the lowest simulated temperature and 



at U=16.5 (closest to the quantum critical point) the 
measured autocorrelation time for kinetic energy at the 
lowest frequency (the slowest mode) was about 0.5 sec 
(see Figure [t]). With the typical simulation time exceed- 
ing three weeks on a single CPU this translates to > 3 
million uncorrelated samples. Note also, that our space- 
time volume in this case contains about 200 correlation 
length volumes. 

We consider the response of the system to a small 
uniform modulation of the optical lattice depth, which 
can be described by a perturbation V which is propor- 
tional to the total kinetic energy of the system, K = 
-JY.<ij>^\^j^ namely 

V = 8{t)K, m = ^-^. (4) 

where 5J{t) is proportional to the lattice modulation am- 
plitude. We consider the response to lowest non-trivial 
order when the response is linear in time and quadratic in 
the perturbation amplitude, which is valid for weak per- 
turbations under otherwise adiabatic conditions of the 
system. Experimentally, this quantity is directly propor- 
tional to the energy absorbed by the system when mod- 
ulation at different frequencies is performed for a fixed 
number of cycles [18]. 

In this regime, we are allowed to consider the correla- 
tion function 

X{iujn) = {K{t)K ^ {K), (5) 

at Matsubara frequencies uOn = 27rTn (n > 0), compute 
it numerically and perform an analytic continuation pro- 
cedure to obtain its spectral density S{uj). The subscript 
iujn denotes that the Fourier transform is taken of the 
corresponding correlation function {K{t)K{0)) in imag- 
inary time. In our simulations, the Fourier transform is 
taken on the flight, i.e., we immediately collect statistics 
for different Matsubara frequencies, which can be done 
much faster than for the correlation function in imagi- 
nary time. The quantity x(^'^n) approaches a constant 
for large Matsubara frequencies, namely minus the ki- 
netic energy —{K). Statistics for the kinetic energy are 
collected separately, and the quantity is afterwards sub- 
tracted from the measurements for the correlation func- 
tion. The correlation function at zero frequency has to 
be discarded because it contains no information on the 
Higgs mode, and is orders of magnitude larger than the 
signal at non-zero Matsubara frequencies. We considered 
system sizes L x L = 20 x 20 and inverse temperatures 
P < 10. The correlation length, determined from the 
one-body density matrix, scales as ^c/<^ ^ 5J/E where E 
is the characteristic energy scale of the problem (T at the 
critical point [Hi Uc = 16.7424, particle- hole gap Ami in 
the MI phase, or the Higgs mode frequency in SF) with 
a the lattice spacing. Except in very close vicinity of 
the critical point, |/7 — /7c| < 0.2, our system size is thus 
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much larger than the correlation length. Since measuring 
correlation functions is one to two orders of magnitude 
more difficult than single particle properties, the increase 
in CPU would be too expensive to study larger system 
sizes, while not altering our conclusions. 



THE DATA AT LARGE IMAGINARY 
FREQUENCIES 
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of the correlation function is indistinguishable within the 
error bars from the one at a lower temperature f3 = 2 (see 
Fig.|9|. We use therefore the tail determined at /3 = 2 for 
the P = 10 case as well, hereby substantially reducing the 
required CPU at low temperatures (we measure at most 
the correlation function for a 100 Matsubara frequencies 
at every temperature). 

The next step is applying a Fourier transform to the 
correlation function in Matsubara domain to the imag- 
inary time domain. Because of the slow convergence of 
the functions l/(ico'n)^, we sum up to 500,000 frequencies 
to ensure that frequency truncation has no detectable 
effect. This correlation function in imaginary time ap- 
proaches, in general, a non-zero constant for r (3/2. 
This non-zero constant will be seen as a delta-peak at 
cj = for S{uj). We subtract this delta- functional con- 
tribution from x{'t) by shifting the data resulting in the 
correlation function Xs{^)^ which is shown in Fig. 10 For 
large r, this function can be fitted to ~ l/r"^, in line with 
theoretical predictions [7, 8 for the scaling behavior of 



the spectral function S{uj) 



at ^ 0, though our 



errors bars in the tail are too large to resolve this law 
unambiguously. 



FIG. 8. (Color online). Monte Carlo data for the correlation 
function x(^<^n) as a function of the Matsubara frequencies 
iun- The solid line is a fit of the form f(x) = ^ + which 
is used for iujn > 50. Parameters are J=l,3 = 2, [7=14 
and (n) — 1. 
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FIG. 10. (Color online). The correlation function Xsi^) as 
a function of imaginary time r over the range r G [0,1.5]. 
For large r, the data are compatible with the asymptotic be- 
haviour l/r^, corresponding to S{uj) ^ for small cj, pre- 
dicted by scaling 7 , 8 . Parameters are J = 1, /3 = 10, [/ = 14 
and (n) — 1. 



FIG. 9. (Color online). Monte Carlo data for the correlation 
function ^(icjn) as a function of the Matsubara frequencies 
iujn for two different inverse temperatures, /3 = 2 and /3 = 10. 
Error bars are shown, but are very small. Same parameters 
as in Fig. |8] 

A typical set of Monte Carlo data is shown in Fig. |8] 
The data for large Matsubara frequencies are fitted using 
the form f{x) = ^ + with a, 6 fitting parameters. 
Only even powers of ujn are allowed by causality and 
symmetry. The fit is used for frequencies iuJn > 50. At 
low temperatures, say P = 10, the high frequency tail 



ANALYTIC CONTINUATION 

We are now in a position to perform analytic contin- 
uation. The spectral density S{uj) > 0{uj > 0) is related 
to the correlation function in imaginary time Xs{'t) by 

Xs{t)= I ckoK{u;,T)S{Lo), (6) 

with the kernel 

Ko{uj, t) = exp(-cjr), (7) 
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seeks the best solution reproducing Xs{^) with additional 
self-consistently adjusting quadratic measures favoring 
smooth non- negative functions. The comparison is shown 
in Fig. [TT] We see that the position of the first and the 
second peak are reliable, just as the integrated weight of 
the spectral density. We do not report any data at fre- 
quencies higher than 20 J in this article. Also till what r 
the function Xs{^) can be computed matters enormously. 

COMPARISON WITH REF. [12] 



FIG. 11. (Color online). Spectral density S{(jj) obtained 
by different analytic continuation methods. Parameters are 
J = 1, = 18, = 1, and (n) 1. The gap is A/J 2.8, 
see Ref. \l7\ in agreement with the onset of spectral weight in 
this figure. The maxent parameters [221123] are — 0-50 and 
a = 0.346 with a flat default model. Error bars determined 
by ALPS [23j on the maxent curve for a few data points are 
shown. 



at zero temperature, and 

K{uj^ r) = exp(— cjr) + exp(— — r)), 



(8) 



at finite temperature T = 1/ j3. Given S{uj)^ the compu- 
tation of Xsij) is easy. The inverse problem is equivalent 
to an inverse Laplace transform, which is exponentially 
sensitive to the data. Only when the function Xs{^) is 
known analytically over a finite interval, analytic con- 
tinuation can be done reliably (analytically), but in the 
presence of statistical noise and grid binning errors, the 
problem is ill-defined. 

Practical schemes seek a compromise between 
smoothening the target function and reproducing Xs{^) 
within error bars. Smoothening is necessary to acquire 
a physically acceptable solution. The degree and type of 
smoothening is never independent of the choices made by 
the practitioner or the used software. However, the bet- 
ter the quality of the data in imaginary time, the more 
acceptable solutions are constrained, and hence the bet- 
ter the analytic continuation can be performed. This 
renders a correct error analysis indispensable, which was 
done here using the alea library of the ALPS project [23] . 
On general grounds, one expects that gaps, the position 
of the first peak and the integrated weight are reliable, 
but that high frequency features are beyond control, and 
that those may have a back-action on the width of the 
peaks. 

We applied three analytical continuation methods to 
our data in order to extract what features are seen by 
different methods, and are hence likely to be physi- 
cal. We used maxent [22] in the ALPS implementa- 
tion ^3] , stochastic inference [28l[^ in the implementa- 
tion of Ref. [30] , and a novel heurisitic procedure which 
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FIG. 12. (Color online). Spectral density of the scalar sus- 
ceptibility according to Ref. [12 for different values of U/J 
based on a large N — 2 expansion in d — 2 dimensions. At 
low cj, the behavior of the scalar susceptibility is uo^ . The am- 
plitude was chosen arbitrarily because of the high frequency 
tail ^ l/cj. The theory predicts a sharpening of the peak 
when U is decreased, and saturation of the response in the 
maximum region when U ^ Uc, i.e. no softening is observed 
at criticality. 

The theory of Ref. [12] predicts for the imaginary part 
of the scalar susceptibility Xppi^) (^o be identified with 
our S{uj))^ 
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/7/8J, 

sV2J{l-U/Ucy^^ 



(9) 
(10) 
(11) 



In this large N expansion, N = 2 for < n >= 1 were set. 
This parametrization follows from App. A in Ref. i24l 
(involving a large occupation number Ansatz. Also a 



8 



16.5 



S(co)t 



16 



14 



12 



QMC 



U=10 



r 



2 4 



8 10 CO /J 



(Orr 



Podolsky et al 



Am,(2Uc-U) 



QMC 



u 



FIG. 13. (Color online). Spectral density S{(jo) for low fre- 
quencies obtained by Monte Carlo and analytic continuation 
for different values of U/J. This should be compared with 
Fig. 
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(see also text). 



constant density of states for the phonons is assumed.) 
and was also used in the analysis of the experiment ^ . 



'\ in line with Ref. [7]- 
a consequence of the 



At low frequencies Xpp{^) ^ ^' 
At high frequencies, Xppi^) ^ ^ 
truncation scheme in the theory. This means, in particu- 
lar, that the normalization integral diverges and thus can 
be chosen arbitrarily for comparison with simulation. 

The crossover from c*;^ to uj behavior at low frequencies 
in Eq. ^ occurs at frequencies ~ w? ~ (1 — U/Uc) 
and constitutes strong violation of the universal low- 
frequency scaling law, S{uj) oc A^~'^^^ F{uj/ A) (with 
A ~ (1 — U/UcY and the correlation length exponent 
V = 0.6717), obtained by dimensional analysis [5]. It 
predicts that after crossover at ~ A the response func- 
tion increases as S{ljj) ^ = J^-^^^^ . 

The (unnormalized) response is shown in Fig. 



12 



and 

should be compared with the (normalized) response com- 
puted by quantum Monte Carlo simulations shown in 
Fig. [13] Equation [9] has the following behavior, 

• For 27 < m (this corresponds to U j J > 15), the 
response is peaked at the Higgs mass. The position 
of the maximum is within the peak width the same 



as the one found in Fig. 13 However, the response 
becomes narrower for smaller values of [// J, unlike 
the Monte Carlo prediction in Fig. [Tsj 

• For 27 > m (this corresponds to 15 < U/J < 
Uc/J)^ the response saturates when tuning closer 



FIG. 14. (Color online). Comparison of the peak maximum 
between QMC (black circles for Higgs mode and blue square 
for the broad peak at U/J — 12) and Eq. |9] (red line). Also 
shown is the gap determined in the Mott Insulator from the 
single particle Green function [ITj by mirror reflection around 
the critical point (thin black line). The maximum agrees for 
U/J — 14, which is where QMC predicts the Higgs mode 
to be only barely discernible. For larger values of U the peak 
maximum in Eq.[9]saturates to the value predicted by 7 alone, 
i.e. in the critical region the peak maximum in Eq. [9]is not 
linked to the Higgs mass. 



to the critical point. No further softening can be 
seen, only a broad response, where the maximum 
can no longer be identified with the Higgs mass. 
Note that the QMC behavior in this regime is differ- 



ent (see Fig. 13): further softening and narrowing 



of the low-frequency Higgs resonance is observed 
suggesting that this peak is part of the universal 
scaling scaling function, and a second peak emerges 
at cj/J ^ 6 (the latter is not part of the theory of 
Ref. HI). 



A comparison for the peak maximum is shown in Fig. 14 



DISAPPEARANCE OF THE HIGGS 
RESONANCE IN THE MOTT PHASE 

On the superfluid side of the transition, our data indi- 
cate that the Higgs resonance is seen only in close vicinity 
of the transition point (within 20%). The control simu- 
lation of the MI phase at U/J = 20.5 confirms that this 
20% figure of merit can be used on both sides of the 
transition, see Fig. [15] 
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FIG. 15. Spectral density in the Mott phase at U/J = 20.5. 



SIMULATIONS IN A TRAP 

In the experiment [18], a particle number of 
{N) = 190(36) was reported. Together with 
the requirement that the density is one in the 
trap center, this corresponds to a trap parameter 
Vc{x,y) = 0.0915((x — Xc)^ + — Vc)'^) in units of the 
hopping for U = 14. Here, {xc^Vc) are the lattice coor- 
dinates of the trap center. The reported temperatures 
T/U ~ 0.1 correspond to an entropy S/N = 0.8 in the 
atomic limit. Assuming adiabaticity at all stages in 
experiment, this would mean that (3 J = 0.6 in the trap 
for U = 14, which is a temperature 50% above Tg. Since 
the passage to the atomic limit is not adiabatic, and 
that substantial energy is added during the 20 cycles 
of the lattice modulation (also due to quantum optics 
effects, as was seen in Ref. [25 , the true temperature 
before the lattice modulation is lower, perhaps by 50%. 
The heating has not been addressed experimentally, and 
therefore a more quantitative estimate of temperature 
is not possible. The rule of thumb T/U ~ const, for 
presently realistic temperatures of the scale of the 
bandwidth and for values of U outside the Bogoliubov 
regime, nevertheless holds, in agreement with previous 
experiments [26l l27l . 
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